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Abstract: A bulk viscosity is introduced in the formalism of modified gravity. It is shown 
2 that, on the basis of a natural scaling law for the viscosity, a simple solution can be found 

for quantities such as the Hubble parameter and the energy density. These solutions may 
incorporate a viscosity-induced Big Rip singularity. By introducing a phase transition in the 
cosmic fluid, the future singularity can nevertheless in principle be avoided. 
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1. Introduction 

- - s Modified gravity has become an active branch of modern cosmology, attempting to give a unified 

9 description of the early (inflationary) epoch of the universe and at the same time intending to account for 

10 the accelerated expansion at the later stages. Useful reviews on modified gravity theories can be found 

11 inRefs. [1-3]. 

12 Most treatises on modified gravity, as well as on standard gravity, assume the cosmic fluid to be ideal, 

13 i.e. nonviscous. From a hydrodynamicist's point of view this is somewhat surprising, since there are 

14 several situations in fluid mechanics, even in homogeneous space without boundaries, where the two 

15 viscosity coefficients, the shear coefficient r] and the bulk coefficient, £ come into play. This means a 

16 deviation from thermal equilibrium to the first order. Such a theory means in effect acceptance of the 

17 Eckart 1940 theory [4]. An important property of the Eckart assumption is that the theory becomes 
is noncausal. By taking into account second order deviations from thermal equilibrium, one can obtain 
19 a causal theory respecting special relativity. Pioneering articles on causal fluid mechanics are those of 
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20 Miiller [5], Israel [6], and Israel and Stewart [7]. A recent review can be found in Ref. [8]. Because of 

21 the assumption about spatial isotropy, the shear coefficient is usually omitted. 

22 Our purpose in the following will be to include the bulk viscosity £ in the modified gravity formalism. 

23 We consider the case when ( is satisfying a scaling law, reducing in the Einstein case to a form 

24 proportional to the Hubble parameter. It turns out that this scaling law is quite useful. We survey 

25 first earlier developments along this line, extracting material largely from our earlier papers [9-11]. 

26 Thereafter, as a novel development we investigate how the occurrence of a phase transition can change 

27 the development of the universe, especially in the later stages approaching the future singularity. (It may 

28 here appear natural to relate such a phase transition with the onset of a turbulent state of motion.) It is 

29 shown that such a transition may in principle be enough to prevent the singularity to occur at all. This 

30 part of the paper, covered in Sect. 4, is a generalization of the viscous/turbulent theory for standard 

31 cosmology recently given in Refs. [12] and [13]. 



32 2. Fundamental Formalism 



The action in modified gravity is conventionally written in the general form 

'F(R) 



S = d xy— g 



2k 2 

where k 2 = 8nG, and where £ mat ter is the matter Lagrangian. The equations of motion are 



(1) 



-g^F(R) + R^F'(R) - V„V„F'(R) + g^UF\R) = n 2 T™ ttc \ (2) 



33 where T™ atter is the energy-momentum tensor corresponding to £ m attcr- 

We shall however in the following not consider the general case, but limit ourselves to the special 
form where 

F(R) = f R a , (3) 

34 / and a being constants. This model has been used before, by Abdalla et al. [14] and others. The 

35 case of Einstein gravity corresponds to /o = 1 and a — 1. This choice appears to be natural from a 

36 mathematical viewpoint, and in our case it will play an important role in connection with the scaling 

37 law for the bulk viscosity; cf. Eq. (25) below. Quadratic Lagrangians were considered also earlier, by 

38 Barrow [15,16]. 

We assume the spatially flat FRW metric 

ds 2 = -dt 2 + a 2 (t)dx 2 , (4) 

and put the cosmological constant A = 0. In comoving coordinates, the components of the four-velocity 
£/ M are U° = 1, U % = 0. Introducing the projection tensor = g^ v + U^U V we have for the energy- 
momentum tensor 

Tfj, v = pUyPv + phftv, (5) 

where p is the effective pressure 

p = p- 3H(. (6) 
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The scalar expansion is 9 = 3a/ a = 3H, with H the Hubble parameter. The shear viscosity is here 
omitted. 

The equations of motion following from the above action are 
1 

~T 

The equation of state for the fluid is written as 



-f o9tlu R a + afaR^R*- 1 - a/oV^iT" 1 + af^DR^ 1 = k 2 T™^ . (7) 



p = w p = (7 - l)p. (8) 

If w = — 1 or p = —p the fluid is a vacuum fluid with strange thermodynamical properties such as 
negative entropies (cf., for instance, Ref. [17]). Recent observations indicate that w = — 1.041q 5 > q 
[18,19]. It has been conjectured that w is a function varying with time, perhaps even oscillatory, and that 
w might have been around at redshift z of order unity [20]. The quintessence region — 1 < w < — 1/3, 
and the phantom region region w < — 1, are both of physical interest. Both quintessence and phantom 
fluids imply the inequality p + 3p < 0, thus breaking the strong energy condition. 

We now consider the (00) component of Eq. (7), observing that R 00 = —3a/ a and R = 6(H + 2H 2 ). 
With T r g attcr = p we obtain 

\f R a - 3af {H + H^R*" 1 + 3a{a - l)f HR a ~ 2 R = n 2 p. (9) 

An important property of (9) is that the four-divergence of the LHS is equal to zero, V i/ T / ^, attcr = [21]. 
This is as in Einstein's gravity, meaning that conservation of energy-momentum follows from the field 
equations. The energy conservation equation becomes 

P+(p + p)3H = %H 2 . (10) 

Differentiating (9) with respect to t and inserting p from (10), we get 

^iUR a + 3af [2H - 3j(H + H 2 )^" 1 + 3a{a - l)/ [(3 7 - l)HR + R]R a ~ 2 



+3a{a - l)(a- 2)f R 2 R a ~ 3 = 9k 2 (H. (11) 

Inserting R = 6 (if + 2H 2 ) we see that this equation for H(t) is quite complicated. We shall be interested 
in solutions related to the future singularity, and make therefore the ansatz 

TT 

H — —7, where X = 1 - BH t, (12) 
X 

where H is the Hubble parameter at present time, and B a non dimensional constant. If a future 
singularity is to happen, B must be positive. 

Before closing this section, it is desirable to comment on stability issues for our ansatz (3) for the 
modified Lagrangian. A theory of modified gravity should admit an asymptotically flat, static spherically 
symmetric solution. Now, we expect that the expression (2) for the complete action, with (3) inserted, 
will not be the full solution. It is reasonable to expect that the modified part will contain also other terms 
so that (2) makes up only a part of the complete action. Nevertheless, it is of interest to ask to what extent 
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(3), when take separately, will behave with respect to the stability requirements. In the Solar system, far 
from the sources, it is known that R fa 10~ 61 eV 2 ; it corresponds to one one hydrogen atom per cubic 
centimeter. [Note that 1 eV= 5.068 x 10 4 cm" 1 .] On a planet, R = R b fa 10~ 38 eV 2 , whereas the average 
curvature in the universe is R fa 10~ 66 eV 2 . According to the stability analysis of Elizalde et al. [22], 
the stability condition for matter is 

F"(R b ) > 0, where R b fa 1(T 38 eV 2 . (13) 

In our case, this means merely that the exponent a in the expression (3) has to be greater than one. The 
stability condition on a is quite modest. 



3. Special Cases 

It is now mathematically simplifying, and physically instructive, to focus on special cases. 



3.1. Einstein Action 



As mentioned above, Einstein's gravity corresponds to f = 1 and a — 1. It means that the 
Lagrangian is linear in R. It is natural to consider this case as a reference case before embarking on 
the nonlinear general situation. 

We first have to adopt a definite form for the bulk viscosity. The simplest choice would be to put 
£ = constant. There are however reasons to assume a slightly more complicated form, namely to put £ 
proportional to the Hubble parameter H. This is physically natural, in view of the large fluid velocities 
expected near the future singularity. Such violent conditions should correspond to an increased value of 
C- We shall take ( to be proportional to the scalar expansion, 9 = 3H, 

C = 3r E H, (14) 

r E being the proportionality constant in the Einstein theory. An important property of this particular 
form, shown in Ref. [23], is that if r E is sufficiently large to satisfy the condition 

X = -7 + 3k 2 t e > 0, (15) 

then a Big Rip singularity is encountered after a finite time t. Even if the universe starts out from the 
quintessence region (—1 < w < —1/3) or (0 < 7 < 2/3), the presence of a sufficiently large bulk 
viscosity will drive it into the phantom region (w < —1) and thereafter inevitably into the Big Rip 
singularity. 

From the governing equations we get 



B = \x, #o = y^ 2 Po, (16) 

where p is the present (t = 0) value of the energy density. The time dependent value p E for the energy 
density according to the Einstein theory becomes 

PE = -^2- (17) 
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We ought here to mention that other forms for the bulk viscosity, more complicated than the form (14) 
above, have been suggested. One possibility is that (, in addition to the term proportional to H, contains 
also a term proportional to a /a. See further discussions on this topic in Refs. [24] and [25]. 

3.2. Modified Gravity Action 

Consider now the modified gravity fluid, for which f and a are arbitrary constants. As before, we 
look for solutions satisfying the ansatz (12). It turns out that such solutions exist, if we model the bulk 
viscosity ( a according to the following scaling law [10,26], 

( a = rj 2 ^ 1 = r a (3H) 2a ~\ (18) 

We see that this scaling fits nicely with our results from the preceding subsection: if a = 1, our previous 
form (14) follows. The time-dependent factors in (11) drop out, and we get the following equation 
determining B, 

(B + 2)"- 1 {9(2 - a)7 + 3[a + 37 + a(2a - 3){3j - 1)]B + 6a(a - l)(2a - 1)B 2 } 



18K2 '?] r„. (19, 



/o V2 

This equation is in general complicated. Let us consider ct = 2,7 = 0asa typical example (recall that 
7 = corresponds to a vacuum fluid). Then we obtain from (19) (r Q — > r 2 ), 

B 3 + 2B 2 -^ = 0. (20) 
o/o 

If the LHS is drawn as a function of B, it is seen that there is a local maximum at B = — 4/3 and a local 
negative minimum at B = 0, irrespective of the value of r 2 . For all positive r 2 there is thus one single 
positive root. This root is viscosity-induced, and leads to the Big Rip singularity. When r 2 increases from 
zero, there is a parameter region in which there are three real roots. Assume this region, and introduce 
an angle (j) G [0, 180°] such that 

/ 243« 2 r 2 \ 

cos0 = — 1 — . (21) 



128 f 

Then the actual value of the root can be expressed as 

£ = ~ + ^cos^ + 240°V (22) 

For instance, if we choose = 120°, the positive solution becomes B = 0.3547. According to Eq. (12) 
this gives the following Big Rip time 

1 1 2 - 819 ^ 
tBR ~ ~rJT ~ ~h~- ( 

t> ti Q -no 

We may also note the general relation for B following from the energy conservation equation (10), 
when p p a ,p -»■ p a , C -»■ Ca> 

B = j_! + ^(mfl, (24) 

2a 2a po 
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Here we used 

3#o \ Po 



(a = T a [—) , Pa = ^, (25) 

65 and for simplicity we used the same initial conditions at t = for the modified fluid as for the Einstein 
ee fluid, poo = Poe = Po, and H 0a = H 0E = H . 

67 4. On the Possibility of a Phase Transition in the Late Universe 

In the preceding we have surveyed bulk viscosity-induced generalizations of modified gravity, 
following essentially the earlier treatments in Refs. [10,11,26]. Our intention in the following, as a 
new contribution, will be to discuss the flexibility that the above model possesses with respect to sudden 
changes in the time development (we will refer to it as phase transitions) in the late universe. The main 
point is the different solutions for B in the governing equation (11) that are possible when the scaling 
ansatz (18) is inserted. We obtain the following algebraic equation for B, for definiteness still assuming 
a = 2, 

5 3 + (2 + | 7 )5 2 + | 7 5-|^ = 0. (26) 
4 * o Jo 

68 This equation generalizes (20) to the case of nonvanishing 7. 

Consider the following scenario: the universe starts out from present time t = and follows the 
equations of modified gravity, with a r 2 -induced bulk viscosity corresponding to a positive value of B. 
That means, the universe develops according to 

X' ^ 2 = T2 \x~) ' P2 = x^' (27) 

69 with X = 1 — BH t. The universe thus enfaces a future singularity at large times. Let now, at a fixed 

70 time that we shall call t*, there be a phase transition in the cosmic fluid implying that the effect from 

71 r 2 goes away. It means that the further development of the fluid will be determined by the 7-dependent 

72 roots of Eq. (26) when r 2 = 0. There are three roots: 

1) The first is B = 0. This is the de Sitter case, corresponding to 

H = p 2 = p*, (28) 

73 where and follow from (27) when t = t*. By assuming that I7I 1, which is of main physical 

74 interest, we see that it is easy to determine the remaining two roots. One of them is 

2) B = —2. This means 

-ff* p* 
H = l + 2H,(t-UY P = [l + 2H*(t-Q] 2 ' (29) 

75 The accelerated expansion is accordingly reversed at t = t*, and the density goes smoothly to zero at 

76 large times. 

3) The third root is B = —37/4, which yields 

H= UpJ^j' P= [l + ^H,(t-Qr (30) 
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The sign of 7 is important here. If the equation-of-state parameter w lies in the quintessence region, 
w > — 1 (7 > 0), then the density of the universe will go to zero for large times, like for the case 2) 
above. By contrast, in the phantom region w < — 1 (7 < 0), the universe will actually move towards a 
Big Rip again, although very weakly so. 

Finally, it is of interest to compare the above results with those obtained in ordinary viscous 
cosmology when the universe, similarly as above, is thought to undergo a phase transition at a definite 
time t*. Such an investigation was recently carried out in Ref. [13] (the one-component case treated in 
Sect. VI). Consider the following model: the universe starts from t = as an ordinary viscous fluid with 
a constant bulk viscosity, 

( = constant = ( , (31) 

and develops according to the Friedmann equations. Assume that the universe is in the phantom region, 
7 < 0. It follows that in the initial period < t < t*, 

h = H ° et/tc cm 

i-§i7i#oa^-i)' 



p e 2t//tc 



L 2 

where t c means the "viscosity time", 



P r1 S| 7 |^ t e ( e */*c_i)]2' (33) 
3 



t c = ( ~« 2 Co ) • (34) 

According to these equations the universe develops towards a Big Rip. Now, after t = t* we imagine an 
era for which 7 turb = 1 + u> t urb > and an equation of state of the form 

Pturb = WturbPturb- (35) 

Here the subscript "turb" refers to our association in [13] of the transition at t — into an era dominated 
by turbulence. 
Then, for t > £*, 

H = 3 -Z— r, (36) 

1 + 27turb-n*(£ — t*) 

P * (37) 



[l + §7turb#*(*-^)] 2 ' 

This means a dilution of the density again, at large times. The Big Rip may thus be avoided, as a result 
of a phase transition in the cosmic fluid. We see that in this sense the behavior is similar in the two cases, 
modified or ordinary, gravity. 

It ought to be made clear that we do not at present have a specific model of the phase transition 
suggested at t = t*. Our association with a turbulent state of motion is however quite natural, on 
the basis of the following consideration: In states of violent local motions of the cosmic fluid near 
a future singularity the transition into a turbulent kind of motion seems physically inevitable, as the 
local Reynolds number becomes then very high. That brings the shear viscosity concept back in the 
consideration, now not in a macroscopic but in a microscopic (local) sense. We expect that there is 
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92 established a distribution of eddies over the wave number spectrum. Most likely this distribution can be 

93 taken to be approximately isotropic, implying the existence of an inertial subrange in which the energy 

94 density is E(k) = aK£ 2 ^ 3 k~ 5 / 3 , where olk is the Kolmogorov constant, e the mean energy dissipation, 

95 and k the wave number. Ultimately, when the magnitude of k reaches the inverse Kolmogorov length 

96 \jr\K — (e/V 3 ) 1 / 4 with v the kinematic viscosity, the local Reynolds number becomes of order unity 

97 and heat dissipation occurs. What we have done above, is to denote the post-transition period t > 

98 conceptually as a turbulent region, where the influence from the bulk viscosity has essentially gone away 

99 and where 7 tur b has become positive, without going into further detail as regards the underlying physical 

100 transition process. 

101 5. Conclusions 

102 Starting from the modified gravity action integral (5), we solved the (00) component of the governing 

103 equation (7) inserting the scaling relation in Eq. (12), H = H /X. The bulk viscosity ( was assumed 

104 in the form (18), generalizing the Einstein-case value (14) frequently used in the literature. It turned 

105 out that the mentioned ansatz (18) permitted solutions in the form (25), corresponding to future Big Rip 

106 singularities. This is thus a Big Rip scenario induced by the bulk viscosity. 

107 In Sect. 4 it was discussed how the future singularity can nevertheless in principle be avoided, if one 

108 allows for a future phase transition in the cosmic fluid where the influence from viscosity goes to zero. 

109 Modified, or conventional, cosmology behave in this sense essentially in the same way. 

no Finally, one may ask to what extent the above theory can be generalized to more complicated forms of 

1 1 1 the modified Lagrangian term than the simple power-law form given in Eq. (3). Such a general expression 

112 for F(R) would then have to be inserted into the field equation (2), and thereafter to be combined with 

113 the energy conservation equation (10). The general case seems difficult to handle, but there may be 

114 special cases that are mathematically tractable and of physical interest. Such an investigation is however 

115 outside the scope of the present paper. At present, the scaling laws like Eq. (25) appears to be closely 
no linked to our basic ansatz (3). 
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